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TITLE OF THE INVErfTlON 

Automatic Multf-Dimenslonal Intravascular Ultrasound Image 
Segmentation Method 

FIELD OF THE INVENTION . 

[OO9I] The present invention generally relates to image segmentation. 

Mo/e specifically but not exclusively, the present invention Is concerned witli an 
• intravascular ulti-asound image segmentation technique* for characterizing blood 

vessel.vascular layers from intravascular ultrasound Image sequences. 

■i 

BACKGROUKD OF THE INVENTION 

[0002] Over the past few years, intravascular ultrasound (IVUS) 

technology has become very useful for studying atherosclerotic disease. IVUS is a 
medical imaging technique that produces cross-sectional images as a catheter is 
pulled-back inside a blood vessel. These images show the lumen but also the 
layered structure of the vascular wall. IVUS provides quantitative assessment of 
the vascular wall, infomnation about the nature of atherosclerotic lesions as well as 
the plaque shape and size such that in clinic, IVUS was rapidly recognized as a 
valuable tool in diagnosis and iii pfe-intervention analysis of atherosclerosis. 

[0003] " The ability to characterize the vascular wall was initially 

demonstrated In 1989. by Gussenhoven et al., in "Arterial wall characteristics 
determined by intravascular ultrasound imaging: An in \^tro study" {J. Am. Coll. 
Cardiol., vol. 14, no. 4, pp. 947-952, 1989). Also, studies of the mld-90s by Mintz 
et al.. in "Atherosclerosis In anglographically 'normal' coronary artery reference 
segments: An intravascular ultrasound study with clinical correlations" {J. Am. Coll. 
Cardiol., vol. 25, no. 7, pp. 1479-1485, 1995). showed, based on IVUS, that 40% 
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of angiographically normal vessel were in fact atherosclerotic. 

[0004] By using IVUS. it was also demonstrated by Colomtx} et a!., in 

"Intracoronary stenting without anticoagulation accomplished with intravascular 
ultrasound guidance' {Circulation, vol. 91, pp. 1676-1688, 1995) tiiat conventional 
stent implantation resulted In incomplete apposition and expansion causing 
tiirombosis, vt^ich had tiie result of dianging the dlnical practice. 

[00051 IVUS Is also expected to play an . Important role in 

atherosclerosis research. For example, as demonstrated' by Nissen et al., in. 
"Application of intravascular ultrasound to characterize coronary artery disease 
and assess the progression or regression of atherosclerosis" (Am. J. Cardiol., vol. 
89, pp. 24B-31B, 2002), IVUS helps to achieve predse evaluation of the disease 
in new progression^regression therapies. Experts also generally agree that IVUS 
imaging adds predous complementary information to angiography whidi only 
shows a projection of the lumen, as taught by Nissen et al., in "Intravascular 
ultrasound : Novel pathophysiological Insights and cuh-ent clinical applications" 
. {Circulation, vol. 103, pp. 604-616, 2001 ). 

[0006] Over the last few years, new signal processing strategies were 

applied to IVUS signals for extracting infonmation on the elastic properties of the 
vascular wall. For sample, a new imaging technique named "intravascular or 
endovascuiar .utoasound elastography" was proposed by de Korte et al., in 
"IntrBvascuiar elasticity Imaging using ultrasound - Feasibility studies in phantoms" 
{Ultrasound Med. Bloh, vol. 23, pp. 735-746. 1997). Recently, Bmsseau et al. In 
"Fully automatic luminal contour segmentation in intracoronary ultrasound imaging 
- A statistical approach" {IEEE Trans. Med. Imag., vol. 23. pp. 554-566. 2004) 
suggested to use a pre-segmentaflon of ttie stiructures of the vascular wall 
identified from IVUS images to help Implementing elastography algorithms. This 
constitutes anottier Important domain of application of IVUS multi-dimensional 
image segmentati'on. 
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[0007] The tomographic nature of IVUS makes 3D reconstruction of the 

vessel wall possible. FurthemiQre, 2D and 3D quantitative measurements of 
atherosqierotic disease such as plaque volume, intima-media thickness, vascular 
remodeling, and lumen area stenosis can be retrieved from IVUS data as 
disclosed by MIntz et a!., In "American college of cardiology clinical expert 
consensus document on standards for acquisition, measurement and reporting of 
intravascular ultrasound studies (IVUS)" (J. Am. QolL Cardiol., vol. 37. no. 5, pp. 
1478-1492. 2001). 

[0008] However, a typical IVUS acquisition generally contains several 

hundreds of images, which has the effect of making analysis of the data a long and 
fastidious task that is further subject to an important variability between intra- 
observers and inter-observers when non-automatic methods are used. These 
aspects generate important constraints against the clinical use of IVUS. Other 
constraints related to the use of IVUS include poor quality image due to speckle 
noise, imaging artifacts, and shadowing of parts of the vessel wall by calcifications. 

[0009] So far, a number of segmentation techniques have been 

developed for IVUS data analysis and were introduced to overcome the 
hereinabove discussed drawbacks. Generally speaking, a portion of these 
techniques are based on local properties of Image pixels, namely with the gradient- 
based active surfaces as introduced by Klingensmith et al., In "Evaluation of three- 
dimensional segmentation algorithms for the identification of luminal and medial- 
adventitial borders in intravascular ultrasound images" (IEEE Trans. Med. Imag., 
vol.. 19, no. 10, pp. 996-1011, 2000) and the pixel Intensity combined to gradient 
active contours as introduced by Kovalskl et al., in.ThreeHdimensional automatic 
quantitative analysis of intravascular ultrasound images' {Ultrasound Med. BloL, 
vol. 26, no. 4, pp. 527-537, 2000). 

[0010] Graph search was also investigated using local pixel features, 

for instance, with Sobel-like edge operator as disclosed by Zhang et al., in "Tissue 
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characterization in intravascular ultrasound images" {IEEE Trans. Med. Imag., vol. 
17; no, 6, pp. 889-899, 1998) and with gradient assoaated to line patterns 
correlation as demonstrated by Von BIrgelen et al., in "Morphometric analysis in 
three-dimensional intracoronary ultrasound : An in vitro and In. vivo study using a 
novel system for the contour detection of lumen and plaque" {Am. Heart J., vol. 
132, no. 2. pp. 516-527, 1996). 

[001 IJ The other portion of the IVUS segmentation work was based on 

more global or region infonnation. For Instance, texture-based morphological 
processing was considered as disclosed by Mojsilovic et al., in "Automatic 
segmentation of Intravascular ultrasound images: A texture-based approach" {Ann. 
Biomed. Eng., vol. 25. no. 6, pp. 105&-1071, 1997). Gray level variances were 
then used for the optimization of . a maximum a . posteriori (MAP) estimator 
modeling ultrasound speckle and contour geometry as demonstrated by Haas et 
al., in "Segmentation of 3D Tntravascular ultrasonic images b^sed on a random 
field model" {Ultrasound Med. Biol, vol. 26, no. 2, pp. 297--306, 2000). 

[0012] In addition, some studies defining only the lumen boundary and 

not using the full IVUS potential can be found in the literature. Still, in 2001, the 
clinical expert consensus froni the American College of Cardiology in the 
hereinabove cited document by IWIhtz et al. reported that no IVUS edge detection 
method had found widespread acceptance by clinicians. 

[0013] Recently, graph search was revisited using other edge filters, as 

disclosed by Koning et al., in "Advanced contour detection for threeklimensional 
Intracoronary ultrasound: A validation - In vitro and in vivo" {Int J. Cardiac Imag., 
vol. 18. pp. 235-248, 2002), 

[0014] Other recent models and methods were proposed, such as 

elliptical template fitting as demonstrated by Weichert 6t al., in "Virtual 3D IVUS 
model for intravascular brachytherapy planning: 3D segnientatlon, reconstruction, 
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arid visualization of coronary artery architecture and orientation" {Med. Phys.,"vo\. 
30, no. 9, pp. 2530-2536, 2003) and multlagent segmentation by Boveniomp et 
ai., in "l\/luitiagent IVUS image Interpretation" (SPIE Proceedings : Medical Imaging 
2003 : Image Processing, vol. 5032, San-Diego, USA. 2003. pp. 619-630). 
However, these new models were again using local pixel or edge infomnatlon and 
they were not taking advantage of the statistical Information of IVUS data (^eckle 
texture). 

[0015] Since image pixels in IVUS have pixel gray values generally 

distributed according to a Rayieigh probability density func^on (PDF) in B-mode 
(brightness modulation) imaging of uniform scattering tissues, as demonstrated by 
Wagner et al,, in 'Statistics of speckle In ultrasound B-scans" {IEEE Trans. Son. 
Unrason., vol. 30, no. 3, pp. 156-163, 1983), it is believed that PDF features can 
be of value for IVUS segmentation. This infomiation is hypothetically more suitable 
for IVUS image analysis, especially when contrast is low between layers of the 
vascular wall. In addition, because the IVUS radio-frequency (RF) mode generally 
provides a better spatial resolution than B-mode imaging, it Is also expected that 
the Gaussian PDF of RF images can be exploited for image segmentation. 

[00161 Since the atherosderolic plaque stnjcture on the vascular wall 

can have an inegular and complex shape that Is rarely elliptical, a fast marching 
method as disclosed by Sethian In 'Level Set Methods and Fast Marching 
Methods: Evolving Interlaces in Compuktiional Geometry, Fluids Mechanics, 
Computer Vision and Materials Science" (2nd ed. Cambridge, UK: Cambridge 
University press, 1999) and by Osher et al., In "Fronts propagating with curvature- 
dependent speed: Algorithms based on hamilton-jacobi fonriulations" (J. CompuL 
Phys., vol. 79, pp. 12-49, 1988), can be used to handle topological changes and 
contour io-egularities generated by IVUS Images. Further, the fact that a fast 
marching method propagates Interfaces in the direction of the boundaries through 
an exhaustive analysis of the propagation- region has the effect, of decreasing the 
variability of segmentation results. 
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SUMMARY OF THE INVENTION 

[0017] More specifically, In. accordancse with the present Invention, 

there Is provided an Image segmentation method for estimating boundaries of 
layers In a multi-layer body, the method Including providing image data of the 
multi-layer body, the Image data representing a plurality of ^Image elements. The 
method further Includes detemnining a plurality of Initial Interfaces coTOspondIng to ' 
regions of the Image datd to segment, and concun^ntly propagating the initial 
Interfeces con-espondlng to the regions to segment and thereby estimating the 
boundaries of the layers of the multi-layer body. Propagating the initial Interfaces . 
Including using a fast marching model based on a probability function describing at 
least one characteristic of the image elements. 

[0018] There Is furthemiore provided an image segmentation method 

for estimating boundaries of "layers In a multWayer body, the method Including 
providing image data of the multi-layer body, the image data representing a 
plurality of image elements. The method further includes determining .a plurality of 
Initial Interfaces corresponding to regions of the image data to segment, and 
concurrently propagating the Initial Interfaces corresponding to the regions to 
segment the regions arid estimate the boundaries of the layers of the multi-layer 
body. Propagating the initial interi^aces Includes using a fast marching model based 
on a gradient function describing at least one characteristic of the image elements. 

[0019] There is furthemiore provided an Image segmentation method 

for estimating boundaries of layers In a pulsating multi-layer blood vessel, the 
method Ihcliiding: providing IVUS Image data of the pulsating multi-layer blood 
vessel, determining initial Interfaces corresponding to the regions of the IVUS 
Image data to segment, dividing wall pulsations of the IVUS Image data.lnto a 
discrete number of phases with adjustable pulsation phase labels, assigning the 
pulsation phase labels to 2D IVUS frames of the IVUS Image data, dividing the 
IVUS image data according to the phases and propagating the Initial interfaces 
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according to a fast marching model by simultaneously estimating a mixture of 
• probability density functions in the IVUS image data for each of the regions to 
, segment and according to each of the phases. 

[0020] There is furthenmore provided an image segmentation method 

for estimating boundaries of layers in a multi-layer body, the method Including: 
providing image data of the multi-layer body', the image data representing a 
plurality of image elements. The method further includes determining initial 
interfaces conresponding to the regions of the image data to segment and 
propagating the initial interfaces according to a fast marching model. Propagating 
the initial interfaces includes, for each region to segment, simultaneously 
estimating a speed function for the propagation of the initial interfaces based on a 
probability function describing at least one characteristic of the Image elements, 
and mapping a fime function of the propagating initial interfaces. 

|]0021] There is furthermore provided a data acquisition system for 

segmenting images by estimating boundaries of layers in a multi-layer body, 
including: a catheter Including a transducer for providing image data representing a 
plurality or Image elements and a data acquisition tool Including: a digitizer in 
communication with the transducer for digitizing the image data, a memory for 
receiving and storing the digitized image data, a calculator for estimating, for each 
of the layers, probability functions describing at least one characteristic of the 
image elements, a processor in communication vwth the memory and the calculator 
for simultaneously estimating the boundaries of the layers of the digitized Image 
data by using a fast marching model based on the estimated probability functions. 

[0022] The foregoing and other objects, advantages and features of the 

present invention will become more apparent upon reading of the following non- 
restrictive description of illustrative embodiments thereof, given by way of example 
only with reference to the accompanying drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS . 
[0023] In the appended drawings: 

[0024] Figure 1 Is a 2D IVUS frame view representing the various 

layers of a blood vessel In Image data which is used In a segmentation method for 
detecting layer boundaries according to an Illustrative embodiment of the present 
invention: 



[00251 Rgure 2 Is another 2D IVUS frame view representing one 

in-egularly shaped layer of the blood vessel shown in Figure 1 ; 

[0026] Rgure 3 is a chart view showng the gray levels present in IVUS 

Image data and the mixture of the probability density functions detected per layer; 

[0027] Rgure 4 Is a longitudinal view generated from IVUS Image data 

showing an operation of a method for detecting the layer boundaries; 

[0028] Rgure 5 Is a 2D IVUS frame view intersecting the longitudinal 

view of Figure 4 at point A; 

[0029] Figure 6 is a flowchart schematically representing simulated 

IVUS image data generated from a plurality of 2D IVUS frames as the one shown 
in Figure 1; 

[0030] Figure 7a Is a cross-sectlonal view of a simulated blood vessel 

reconstructed according to a segmented 2D IVUS frame; • 

[0031] Figure 7b Is a 2D IVUS frame view of the simulated blood vessel 

of Figure 7a and generated by the method shown in Rgure 6; 



wo 2005/04S190 



9 



PCT/CA20<My001970 



[0032] Figure 7c is a 2D IVUS frame view showing the segmentation 

results on the 2D iVUS frame view of Figure 7b, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 

[0033] Figure 7d Is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 7b, with the boundaries detected by 
the 3D gradient fast marching method; 

[0034] Figure 8a is another cross-sectional view of a simulated blood 

vessel reconstructed according to a segmented 2D IVUS frame; 

[0035] Rgure 8b is a 2D IVUS frame view of the simulated blood vessel 

of Figure 8a and generated by the method shown in Figured; 

[0036] Figure 8c Is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 8b, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 

[0037] Figure 8d is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 8b, with the boundaries detected by 
the 3.D gradient fast marching method; 

[0038] Rgure 9a Is a 2D IVUS frame view similar to the one shown in 

Figure 1 ; 

[0039]. Figure 9b ls a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 9a, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 



[0040] 



Figure 9c Is a 2D IVUS frame view/ showing the segmentation 
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results on the 2D IVUS frame view of Rgure 9a, vwth the boundaries detected by 
the 3D gradient fast nriarching method; 

[0041] Figure 10a is a 2D IVUS frame view similar to the one shown in 

Rgure 1; 

[0042] Rgure 10b is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Rgure 10a, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 

[0043] Rgure 10c is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Rgure 10a, with the boundaries detected by 
the 3D gradient fast marching method: 

[0044] Rgure 11 Is a longitudinal view showing a volumic 

reconstruction of the vessel layers detected according to the fast-marching method 
based on probability density functions; 

[0045] Figure 12a is a 20 IVUS frame view representing the various 

layers of a blood vessel on simulated RF Image data; 

[0046] Figure 12b is a 2D IVUS frame view showing the segmentation 

results on the RF image data of Rgure 12a, with the boundaries detected by the 
3D fast-marching method based on probability density functions; 

[0047] Figure 1 3a is. a detailed schematic view showing a first example 

of propagation area for detecting a layer boundary using the fast marching method: 

[0048] Figure 13b Is a detailed schematic view showing another 

example of propagation area for detecting a layer boundary using the fast 
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marching method; 

[0049] Figure 14a Is a 2D IVUS frame view which is undersampled vwth 

respect to a typical 2D IVUS frame and which may be used according to a second 
. illustrative embodiment of the present invention; " 

[OdSO] Figure 14b is another 2D IVUS frame view which is 

undersampled with respect to a typical 2D IVUS frame but with a higher resolution 
than the 2D IVUS view shown in Figure 14a; 

[0051] Figure 14c is another 2D IVUS frame view which is 

undersampled with respect to a typical 2D IVUS frame but with a higher resolution 
than the 2D IVUS view shown in Rgure 14b; 

[0052] Figure 14d is a typical 2D IVUS frame view with a higher 

resolution than the 2D IVUS views shown in Figure 14a to 14c; 

[0053] Figure 1 5 is a detailed schematic view showing template regions 

of a vessel for detecting layer boundaries in a method according to a third 
Illustrative embodiment of the invention; 

[0054] Figure 16 is a cross-sectional longitudinal view of an IVUS 

image data showing a sawtooth artifact typically caused by pulsations of blood 
vessels; 

[OOSq Figure 17 is a chart view showing layer areas resulting from 

segmentation detected by a 3D fast marching method based on probability density 
functions according to a fifth illustrative embodiment of the present invention; 



Figure 18 Is a flowchart representing a method for detemnlnlng 
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the initialization contours according to tiie .first illustrative embodiment of the 
present invention; 

[0057] Figure 19 Is a block, diagram representing an In-vjvo IVUS data 

acquisition and processing device according to the first Illustrative embodiment of 
the present Invention; 

[0058] Figure 20 Is a flowchart representing an IVUS method for 

initializing propagating interfaces generated from lower resolution segmentation 
results according to the second illustrative embodiment of the present invention; 

[0059]- Figure 21 is a flowchart representing an exploration method 

from a wide propagating area at low resolution to a reduced propagating area at 
high resolution according to the second illustrative embodiment, of the present 
invention; 

[0060] Figure 22 is a flowchart representing a template region 

searching method for automatically finding Initial Interfaces of the layers according 
to the third Illustrative embodiment of the present invention; 

[0061] Figure 23 Is a flowchart representing an automatic estimation 

method of the probability density function mixture parameters based on the 
iterative conditional estimation acpording to the fourth illustrative embodiment of, 
the present invention; and 

[0062] Figure 24 is a flowchart representing the method for using layer 

pulsation assessment In the boundary detection segmentation process of the fast- 
marching method according to the fifth illustrative embodiment of the present 
invention. 
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DETAILED DESCRIPTION OF THE ILLUSTRATIVE EMBODIMENTS 

[0063] The non-restriGllve illustrative embodiments of the present 

invention relate to a method and device for concurrently estimating boundaries 
between the plurality of layers of a blood vessel from l\AJS Image data. The 
method and device Involve a segmentation of the IVUS image data by propagating 
interfaces in each layer to be estimated from initial Interfaces that are generated 
from the IVUS Image data. The technique for estimating the boundaries of the 
various layers uses a fast marching method based on probability functions, such 
as for example a probability density function (PDF) or gradient function to estimate 
the distribution color map of images, such as for example to estimate the gray 
levels or the muIti-<:olored levels in images, 

[0064] The following description is organized as follows. First of ail, a 

PDF estimation technique for the different vessel layers will be presented. Then, 
an IVUS 3D fast marching method based on the estimated PDFs and based on the 
gray level gradient will be explained and followed by an initializing technique. 
Finally, segmentation results on experimental B-mode data, simulated B-mode and 
simulated RF data and will be reported and discussed. 

[0065] IVUS images are generally provided from an ultrasound 

transducer at the tip of a catheter that is pulled bade Inside a blood vessel and 
produces a plurality of IVUS 2D frames. A typical IVUS 2D frame Is illustrated in 
Figure 1. As illustrated, the 2D frame of Figure 1 shows the catheter and some 
layers of the biood vessel such as, for example, the lumen, the Intima and plaque, 
the media and the sun-ounding tissues. Figure 2 illustrates how the boundary of the 
. lumen may be In^egulariy shaped. 

[0066] The IVUS 2D frames are ultrasonic images made from a 

plurality of pixels generally colored with various shades of gray. In B-mode 
(brightness modulation) or RF-mode (radio-frequency) imaging, such as for 
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example in IVUS data, a Rayleigh or a Gaussian probability density function (PDF) 
can be used, respectively, to mode) the color rpap distribution of the ultrasonic 
speckle pattern In a tinlfomi scattering tissue. When more than one layer of tissue 
is present, the color map distribution of a whole IVUS Image data can then be 
modeled by a mixture of Rayleigh or Gaussian PDFs, depending on the mode 
selected on ttie instrument. 

[0067] The illustrative embodiment that follows generally considers 

IVUS Bnnode imaging, but one ordinary skilled in the art will easily understand that 
similar equations can be prowded for Gaussian PDFs if the RF-mode is 
considered. For more details, see Hasfie et al.. In The elements of statistical 
learning. Date, mining, inference and prediction" (New York, USA: Splnger. pp. 
236-242,2001). 

[0068] In this first illustrative embodiment, a Rayleigh probability 

density function (PDF) Px{x) models the gray level color map distribution by using 

a parameter a^, where x Is the gray level taking values situated, for example, in 
the range [1. .... 256]. In this particular example, the Rayleigh probability density 
function (PDF) is given by equation 1 : 



[0069] ■ ' ' a' \, 2a ; ' t"^) . 

[0070] with *»«^ > 0^ and the variance =a^i^-fc)/2 , 

[0071] IVUS data are modeled by a mixture of M Rayleigh PDFs 

(con-esponding to iVI different layers of the blood vessel) with parameters 

0 = {(ffl^.aj)}", -^here Is the proportion of the 7* component of the mixture of 
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the M Rayleigh PDFs, so that The global data PDF mixture then 

becomes : 



[0072] p^^ix 



= (2) 



[0073] To describe the PDF mixture for the globaMVUS data; the 

parameters of each PDF composing the mixture need to be estimated. In 

IVUS data, the occurring probability of the gray level values x, or observed data, 
can generally be measured by computing the image histogram, as shown in Figure 
3, but the blood vessel layer to which each pixel of an IVUS Image belongs is 
generally unknown or hidden for images that are not already segmented. 

[0074] The Expectation-Maximization algorittim (EM) Is an Iterative 

computation technique of maximum likelihood estimates for incomplete data, as 
presented by Dempster et al., in "Maximum likelihood from Incomplete data via ttie 
EM algorithm" (J. Roy. Stat Soc. B, vol. 39, no. 1, pp. 1-38, 1977), which can be 
used to provide tiie unknown parameters or hidden information of the pro|)ability 
density functions (PDFs). Because ttie IVUS data are Incomplete in tenns of 
maximum likelihood estimation, the EM algorithm can be applied to evaluate Oie 
missing or hidden mixture parameters of tiie Rayleigh or Gaussian PDFs. 

[0075] The EM algoritiim therefore helps to describe the global data 

PDF mixture because 0, a mixture parameter . maximizing tiie likelihood of 
p(3\&) , cannot be solved analytically. A hidden variable Y representative of the 
tissue class (vascular layer to which flie pixel belongs) and taking values situated 

in ttie range [1 M\, must be Introduced at tills point. The log-likelihood of the 

joint distribution of (x,y)={(;c,,:)^,)£„ where N represents tiie size of ttie IVUS 
data, is: , • 
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[0076] lodpx.Y^ix,y\e))=f^logp{yM^i\yi^®) (3) 

[0077] The first step of the EM algorithm is called the Expectation Step 

which calculates the cost function .j2(e,0')=-^Y&og(P(x,Y|©))lX,©'], the 
expected value of the log^ikelihood of (X,Y), the joint distribution, given the 
observed data X and ®' = {(<o^ a previous estimate of the PDF mixture 
parameters. 



[0078] Tiie next step is 'to determine a new estimate © of the PDF 

mixture parameters by maximizing Q{@,&) with respect to parameters ©•; this 
operation can now be perfonfned analytically. 

[0079] The detailed PDF parameter estimation procedure via the EM 

algorithm is therefore: 

[0080] • initialize ©' , the previous estimate of PDF mixture parameters. . 

[0081] • Expectation Step: 

[0082] Evaluate the cost function : 

[0083] e(0,©')=A[log(P(X,Y|®))lX,0'] (4) 



[0084] 



>3 M 



(5) 
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[008^ 



Calculate /jCy/ = / I Xi,©') = u 



according to Bayes rule, and using the. previous parameter estimate 0' and 
Equation 1. 



[0086] 



Maximization Step: 



[0087] 

parameters 



Calculate (§), the new estimate of tiie PDF mixture 



[0088] 



( ( ^ \\ 

fi(©.0')+^ l-E'^y 

I, • \ . M J) 



[0089] 



N ir. 



(6) 



[0090] 
tol. 



where A = N is a Lagrangian making the sum of the ©y equal 



[0091] fly* = argmax^ fi(©,©') 



[0092] 



a, = 



(7) 



[0093] • Jf ©=»^®', update previous estimate ©' = ©, and repeat the 

Expectation and Maximization steps. 
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[0094] In summary, the EM. algorithm maximizes the likelihood of the 

joint distribution of the observed and hidden data by estimating the posterior 
distribution with Pypi,&(y\x,®'). An interesting property of tiie EM algorithm is that 
it is guaranteed tiiat the likelihood of the observed data X increases at each 
iteration. 

[0095] For computation effidency, the EM algorittim is generally 

applied to a randomly drawn subset of tiie observed data X. which are, Iri this 
case, a portion of the pixels fnDm the whole iVUS data. For Instance, the subset 
size may be about 400 000 pixels when a complete IVUS pullback generally 
contains over 80 000 000 pixels. 

[0096] EM algoritiims are otherwise well known to those of ordinary 

skill In the art and. accordingly, will not be further described in the present 
specification. 

[0097]. ' The estimated gray level PDFs of the blood vessel layers can 

then be used to establish a segmentation model in the fast marching framewori<. 
The fast marching method is derived from the level-set model disclosed by SeHiian 
in 'Level Set Methods and Fast Marching Methods: Evolving Interfaces in 
Computational Geometry. Fluids Mechanics, Computer Vision end Materials 
Science' (2nd ed. Cambridge. UK: Cambridge University press. 1999) and by 
Osher et al., in "Fronts propagating with curvatur&KJependent speed: Algoritiims 
based on hamitton-Jacobi formulations" (J. Comput. Phys., vol. 79, pp. 12-49. 
1988). The fast marching metiiod helps to follow interface propagation. 

[0098] In ttie level-set model approach, an Initial interiiace is defined as 

the zero level of a function of a higher dimension than the Interface. The value 
i^{x) of a point X = (x1. x2. .... x„) e 91" is the distance between that point and the 
initial Interface. The function 4> moves in its normal direction according to a speed 
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function F. The evolution of function ^ is given by the following Equation 8 with 
initial interface ^(x,^ = 0). 



[0099] ■^+i?|V«5| = 0. (8) 



[0100] The level-set model Is applicable to Image segmentation by 

interpretihg image boundaries as the final position of the propagating interface, as 
disclosed by Sethlan in "Level Set Methods and Fast Marching Methods: Evolving 
■Interfaces in Computational Geometry, fluids Mechanics, Computer Vision and 
Materials Science" (2nd ed. Cambridge, UK: Cambridge University press, 1999) 
and by R. MalladI et al., In "Shape modeling with front propagation:. A level set 
approach" (IEEE Trans. Pattern Anal. Machine Intell,, vol. 17, no. 2, pp. 158-175, 
1995). 

[0101] To achieve this, the speed function F is defined In tenns of 

image or shape features and should become close to zero when the propagating 
interi'ace meets with image boundaries. Since the speed value is near zero, the 
propagating Interface stops on the image boundary, which generally ends the 
segmentation process, 

[0102] Fast marching is a particular case of the level-set model. It , 

consists of an interi'ace evolving under a unidirectional speed function. In this case; 
for a positive speed function, the propagating interface must be inside the region to 
segment (or outside for a negative speed function), because the propagating 
Interface does not explore the region located inside the initial interface. 

[0103] In the fast marching formulation, the evolution of the 

propagating Interface is expressed in temns of th? arrival time T(x) of the Interface 
at point (x). The function T satisfies the following Equation 9, stating that the arrival 
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time difference between two adjacent pixels increases as the velocity of the 
interface decreases. 

[0104] \77\F=i. (9) 

I 

[0105] The propagation of the interface is done via the construction of 

the anival time function T(x). The construction aigorithm. as disclosed by Sethian 
In "A fast inarching level set method for monotonlcally advancing fronts.' 
.{Proceedings of the National Academy of Sciences of the United States of 
America, vol. 93, pp. 1591-1595. 1996), selects the interface point x having the 
smallest arrival time and calculates the arrival times of its neighbors. This is 
repeated until the interface has propagated across the whole image or until the 
propagating Interface Is considered stationary (when the time gradient is 
sufficiently high). 

[0106] The level-set and fast marching equations are Independent of 

the interface dimension. On a discrete 3D grid, neighbors' anival times are 
updated by solving the following approximation of Equation 9: 

[0107] ■~-mtac(p;^,T-D:jj,T,6f 
[01081 +msK(D;j^T,-D:jj;rfiJ 



[0109] 



^■xB2x{p;^j,-D;!,M (10) 



[01 10] For the x dimension, 
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10111] i)*;^r=±(r^j^-7;^^)/A 

[0112] where A is the grid element size and {i,j,k) is the 3D position 

of the point having its arrivai time calculated. Similar definitions apply for 
ZJ^jjr and 2>5?t3r , in the y and z dimensions. 

[0113] As stated hereinabove, since multiple contours (lumen, intima 

and media) must be identified on the IVUS data, image segmentation is 
simultaneously done via a multiple interface extension of the fast marching method 
as disclosed by Sifakis et al., in "Bayesian level sets for image segmentation" (J. 
Vis. Commun. Image R.,' vol. 13, pp. 44-64, 2002). A speed function is then 
defined for each propagating interface and the T map is built by selecting the point 
with the smallest amval time value from ail propagating Interfaces. 

[0114] Therefore, the fast marching method with multiple propagating 

interfaces enables simultaneous segmentation ;of different layers of the blood 
vessel. The multiple interfaces directly depict the layered «tmcture of the blood 
vessel and provide that the boundaries do not overlap. 

[0115] In the PDF-based fest marching method, each interface 
associated to a vessel layer evolves at a velocity defined in temis of the PDF 
of the corresponding anatomical structure. The propagation speed of the Interface 
wei, where i Is the set 1, 2 of the JV^ evolving interfaces, Is given by 

Equation 11. 



[0116] F„{hhk)= 



1 rj,. log^ 



(11) 
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[0117] /, is the gray level value of pixel s at position {hj\k) in image 

7 , P„{l,) and P,{l,) are the measured occuning probabilities of pixel 7, in region ' 
TO and /, respectively. Because the occumng probability is more significant for a 
region than for a single pixel, the speed function is calculated over a certain 
number of neighbors, which are pixels located around position (1,7,*), such as 
for example, the 26-connected pixels around position {i,j,k). According to 
Equation 11. the m Interface velocities vwll usually be positive and take higher 
values when inside a region having a grayscale distribution dose to P„ . 

[0118] As the propagating interfaces approach the boundaries of the 

blood vessel layers, neighbors start to be distributed under other components of 
the PDFs as stated hereinabove, which has the effect of generally Increasing 
P,{l,) and decreasing P„(7,) and therefore, decreasing the interface speed. The 
velocity function of Equation 11 has a general fonn that may be used with any 
.types of PDF and provides neighboriiood averaging. 

[0119] When used for multiple propagating Interfaces, the fast 

marching segmentation method ends when all adjacent propagating interfaces 
meet with their respective boundaries. Propagating Interfaces thus evolve until the 
anrival time 7 map is completely built. 

[0120] Since the gray level gradient is a widely accepted Image feature, 

comparisons can 'be • made between the hereinabove disclosed PDF 
implementation of the fast marching segmentation and a gray level gradient, 
implementation of the fast marching segmentation. In the latter case, the speed 
function is given by: 
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[0122] where is a Gaussian smoothing filter of standard deviaition 

(T. The speed function of Equation 12 generally propagates interfaces faster on 
low gradient regions. 

[0123] As stated hereinabove, the fast marching segmentation method 

generally requires that the Initial interface is located inside the region to segment. 
This requirement can be achieved with an initialization procedure, in which 
initialization contours are manually traced with respect to data extracted along 
, longitudinal planes of the IVUS data. It can also be automatically perfoniied by 
considering a priori infomnation on the whole IVUS data set, as y/M be further 
described hereinafter. 

[0124] Generally speaking, the step of selecting data along longitudinal 

planes within the IVUS data Is used, instead of using data from a single 2D IVUS 
frame, since longitudinal- planes are able to provide information about the whole 
series of data along the length of the blood vessel. Further, the number of 
manually or automatically traced Initialization contours on the longitudinal plane Is 
Independent of the number of IVUS 2D frames. 

[0125] Initialization contours may be drawn from different numbers of 

longitudinal planes along the blood vessel. As an example, 3 longitudinal planes 
taken at equally spaced angles over 360 degrees may be selected to cut the IVUS 
data volume. The Initialization contours provide reference points for generating the 
set of Initial interfaces on each IVUS 2D frame, for each of the layers to be 
estimated. This is generally accomplished by attributing respecHve reference 
points to the IVUS 2D frame corresponding to each initialization contour points. 

[0126] In the illustrative embodiment of Figures 4 and 5, a reference 

point A is taken along to one initialization contour B estimating the lumen boundary 
on the longitudinal plan© shown in Rgure 4 (Operation 171 of Figure 18). For 
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instance, the reference point A of Figure 4 Is transposed on the IVUS 2D frame of 
Figure 5 to help, generate the initial interface of the lumen on the IVUS 2D frame 
(Operation 172 of Figure 18). More than one reference points generally coming 
frorh more than one longitudinal plane are then transposed on the IVUS 2D frame 
(Operation 173. of Figure 18). This Initialization step is generally done for each 
boundary Jayer of the blood vessel which needs to be estimated. 

[0127] For each IVUS 2D frame-, slightly shrunk splines passing 

through these reference points are computed -and used to generate the Initial 
Interfaces (Operation 174 of Figure 18)! Thereforer using this, procedure, only a 
few longitudinal planes containing IVUS data volume Information are required to 
manually generate the Initialization contours In the longitudinal planes and thereby 
to generate the Initial Interfaces (Operation 175 of Figure 18) required to Initialize 
the concuirent segmentation of multiple vessel layers (Operation 176 of Rgure 
18), over several hundreds of images which Is a typical number for a typical IVUS 
data. . • ■ . 

[01281 The initial longitudinal contours can also altematively be 

restarted from any previous reference points (Operation 177 of Figure 18). In this 
manner, the user can explore, on-line and easily, sections of an IVUS data that 
were more difficult to interpret on longitudinal planes. 

[0129] Experimental Testings 

[0130] Experimental testings of the hereinabove proposed non- 

restrictive Illustrative method of . Rgure 18 were conducted on a total of 8 In-vlvo 
IVUS pullbacks (600 frames/IVUS data) from diseased superficial femoral arteries. 
These experimental testings were performed on 6 different patients before 
undergoing balloon angioplasty. B-mpde IVUS data were acquired with a data- 
acquisition system, such as for example a Jomed equipment (In-vlslon gold, 
Helsingborg. Sweden), iislng a 20 MHz ultrasonic transducer (181 In Figure 19). 
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I 

IVUS 2D. B-mode images of size 10 x 10 mm were digitized (digitizer 182 in Rgure 
19) on 384 x 384 pixel matrices and stored in a memory (183 In Rgure 19) asing 
the DICOIVI standard. The acquisition was done at a 10 images/sec frame rate and 
the catheter puliback velocity was set to 1 mm/sec. generating 0.1 mm thick 2D 
slices. Acquisition parameters were set by a clinician to optimize image quality; 
more spedficaliy, the gain varied from 46 to 54 and the grayscale look-up table 
was set to 5. fmage acquisition was not ECG-gated. 

[0131] To evaluate the robustness of the PDF mixture parameter 

estimation of Oj and a^, the hereinabove described EM algorithm was run 10 

times in a calculator (184 in Rgure 19) for one IVUS catheter puliback, with 
different subsets of pixels at each run of the algorithm. Average PDF mixture 
parameters and standard deviations were calculated for the detected Rayleigh 
PDFs. Since pixel subsets were from the same IVUS data, PDF mixture 
parameters were expected to generally converge toward the same values. 

[0132] * Once this robustness validation was completed, the EM 

algorithm was applied at the beginning of each segmentation, because PDF 
mixture parameters are spedfic to each IVUS data, as' gain and other parameter 
settings are different between each IVUS data, and as echogenicity of each layer 
is variable from one patient to the other. The detected PDF mixtures were 
composed of 4 distributions (lumen, Intima, media, and sunrounding tissues), but a 
skilled reader will easily understand that the EM algorithm fs general and may 
estimate more PDF distributions of heterogeneous vessel layers If required. 

[0133] Testings were conducted on in-vivo blood vessels and 

numerical simulations of blood vessel IVUS data. The in-vivo B-mode IVUS 
images were segmented with 3D multiple interface fast marching using 
automatically detected gray level Rayleigh PDFs and, as a comparison/using the 
gray level gradient. All catheter pullbacks were segmented three times with both 
3D methods using different sets of initial contours. Lumen, intima (plaque), and 
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media boundaries were obtained. To quantify tiie variability of boundary detedion 
under various initializations, average and Haussdorf point-to-point distances, as 
disclosed by Chalana et al., in "A methodology for evaluation of boundary 
detection aigorittims ofi medical images' (IEEE Trans. Med. Imag., vol. 16, no. 5, 
pp. 642-652, . 1997), between resulting contours from different Initial contour sets 
were calculated. Haussdorf distance represents the worst case since it generates 
the maximum distance between different segmentation results. Average and 
Haussdorf distences directly depict polnt-tOi>oint contour variations. 

[0134] Detected boundaries from a wliole IVUS catheter pullback 

represent a blood vessel in 3D that can be reconstructed. Reconstruction of the 
lumen and media contours was made from a sinriple, smoothed contour stack (see 
Figure 11). 

[0135] in addition to the above-described In-vlvo validation of the 

illustrative embodiment of the segmentation method, numerical simulations of 
IVUS data were conducted to evaluate segmentation accuracy. Since the exact 
geometry of the simulated data is generally entirely known, direct perfonnance 
calculations of the detected boundary with respect to the. exact geometry of the 
simulated data can be obtained. The simulated B-mode IVUS data were first 
segmented using the same procedure as for the In-vlvo data, also including 3 
different sets of initial longitudinal view generating Initialization contours. Lumen, 
Intima (plaque), and media boundaries were obtained. Average and Haussdori^ 
point-to-point distances between detected contours and ground truth boundary 
position were calculated for the segmentation results from each set of Initial 
contours. 

[0136] Because the simulation method described In Figure 6 allows 

synthesizing both RF and B-mode IVUS data, the hereinabove described fast- 
matching segmentation mettipd was also tested by using automatically detected 
gray level Gaussian PDFs from the 3D simulated RF images. 
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[0137] The Image-formation mode! that was used to simulate IVUS 

data (echograms) Is detailed by Maurice et al., In 'Adapting the Lagrangian 
spedde model estimator for endovascular eiastographyi Theory and vaHdation with 
simulated radio-frequency data" (J. Aaoust. Soc. Am., vol. 116, pp. 1276-1286, 
2004). Under assumptions sudi as space-lnvarlance of the imaging system, IVUS 
images were modeled by a convolution operation between the point-spread 
function, which is .the equivalent radio-frequency Image of a single ultrasound 
scatterer, and the function describing the acoustic impedance mismatch of each 
scatterer of the simulated tissue structures composing the IVUS data. In other 
words, the point-spread function expresses the Intrinsic characteristics of the 
ultrasound Imaging system. •» 

[0138] The implementation of the image-formation model was made for 

a 20 MHz transducer with a 60% bandwidth at -3 dB and a beam width of 0.1 mm. 
For the purpose of these simulations, the media was selected 2 times more 
echpgenic than the lumen; the plaque, 1.5 times more echogenic than the media; 
and. the sun-ounding tissues 2 times more echogenic than the media. The 
echogenicity can be seen as the image intensity reflecting the acoustic impedance 
mismatch of the scatterers; The signal to noise ratio (SNR) was set to 20 dB. 

[0139] Figure 6 shows the image fomnatlon model used to simulate^ RF 

and B-mode IVUS data. From real 2D in vivo IVUS images, as the one shown in 
box C, the segmented contours or yessel boundaries (lumen, plaque of the intima, 
media) are created (box D) fi-om manually traced contours on the IVUS 2D image 
in box C. Box E shows the function z(x,y) representing ttie acoustic Impedance 
variations within the 2D frame from box D, and box F shows, a function z(r.q>) 
expressing the same acoustic impedance variations mapped wltiiin the 2D frame 
in polar coordinates. Box G shows the polar point spread function h(r, cp) with a 
beam width ttiat increases with deptti and element H is a 2D-convolution operator. 
Processing of function z(r,cp) (box F) wl«i the polar point spread function h(r, q>) 
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(box G) through the 20-convolut)on operator (element H) produces a simulated 
polar radlo-finequency (RF) image l(r,g)) (box J). Box K illustrates the polar B-mode 
Image laCr, qp) that was computed using the Hllbert transfomi (see element L) of 

l(r,<p), as presented by Kallel et al. in "Speckle motion artifact under tissue 
rotation Z {IEEE Trans. Ultrason., Ferroelect, Freq. Contr., vol. 41, pp. 105-122, 
1994). Box M shows the. Cartesian B-mode image or simulated IVUS image iB(x. y) 

calculated from the polar B-mode image jeCr, qp). This simulation strategy was 
repeated for the whole Image data of an IVUS catheter pullback within a diseased 
superfidal femoral artery. 

[0140] Results and Discussions 

[0141] As stated hereinabove, the EM algorithm was applied 10 times 

on 1 • IVUS catheter pullback to evaluate the robustness of the PDF mixture • 
parameter estimation. At each run, PDF parameters were estimated on different 
pixel subsets of the same IVUS data (subsets contained approximately 400 000 
pixels). Average mixture parameters for each detected Rayleigh PDF are shown in 
the following Table 1. An example of automatically detected Rayleigh PDF mixture 
is also shown in Figure 3 In dotted lines for each layer components, with the 
hereinbefore presented IVUS pullback gray level histogram of the whole data set. 

Table I 



Laver 
Comoonent 






Lumen 


29.40±0.10 


0.6!456±0.0021 


Intima and 
Plaque 


20.96±0.50 


347.70±13.02 



wo 2005/048190 PCiyCA2004/001970 

29 



Media 


13.55±0.14 


22.68±0.53 


Surrounding 
Tissues 


36.09±0.67 


2294.58±34.01 



[0142] Table I shows tfiat small variations were found between different 

runs of the EM algorithm. It can be stated that mixture detection of the various 
boundary layers is a robust and stable process with standard deviations of and 
a! ranging from 0^3% to 3.7% for several runs of the EM algorithm applied- on 
different pixel subsets of the IVUS catheter pullback. The EM algorithm Vi^as thus 
applied to the 8 available IVUS catheter pullbacks to study PDF variability between 
different patients. The results are shown in the following Table II. 



Table li 



Laver 
Comoonent 


U)(%) 


a! 


Lumen 


18.82±10.44 


5.52±12.50 


Intima and 
Plaque 


27.81±14.54 


1052.40±1405.97 


Media 


15.87±3.61 


339.46±817.80 


Sun-ounding 
Tissues 


37.50±13.82 


2580.49±664.49 



[0143] Begause of Instrument settings and echogenicities specific to 

the different vascular structures for a given patient, Table .11 emphasizes the 
generally high variability between mixture parameters of distinct IVUS catheter 
pullbacks. These results suggest fliat the EM algorithm Is capable of fitting various 
Raylelgh PDF mixtures from different patients. 
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[0144] The numerically simulated and in-vivo IVUS images can then be 

segmented with 3D multiple interface fast marching methods using automatically 
detected gray level PDFs and gray level gradient for comparison purposes. For the 
experimental testing, all IVUS catheter pullbacks were segmented three times with 
both 3D methods using different sets of initial interfaces obtained from the 
initialization contours generated from the longitudinal planes. 

[0145] The results obtained for the simulated segmentation of the IVUS 

images with the detected gray level Rayleigh PDFs method and with the gray level 
gradient method are shown In Figures 7a to 8d. Figures 7a to 7d are concerned 
with a first blood vessel geometry, shown in Figure 7a, and Figures 8a to 8d are 
concerned with a second blood vessel geometry shown In Figure 8a. 

[0146] Figures 7b and 8b respectively show the simulated IVUS cross- 

sectional 2D B-mode Images of the first and second examples of a simulated blood 
vessel, which can be obtained from the above-described method schematically 
illustrated in Figure 6. In Figures 7c and 8c, the detected boundaries of the lumen, 
the thickened intima and media are estimated with the detected gray level PDFs 
method and, In Rgures 7d and 8d, with the gray level gradient method. 

[0147] The typical simulated IVUS segmentation results shown in 

Rgures 7a to 8d Illustrate that detected boundaries were very close to the blood 
vessel layer boundaries. They also reveal that the external boundary of the media 
is smoother with the PDF fast marching than the gradient-based method, but that 
the lumen, which can have a rougher surface, was detected with sufficient details. 
Gradient methods seemed to trace speckle contours on object boundaries, 
because speckle generally has high gray level Intensity differences. 



[0148] The following Table III includes the results, of the average 

distance (AD) and Haussdorf distance (HD), which Is the maximum distance, 
between the estimated boundary and the true boundary In mm, between detected 
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boundaries obtained from different initialization steps and ground truth values (true 
boundaries) obteined from the simulated geometry: In this table, FMM refers to the 
fast marching segmentation method. Symbol * indicates a statistically significant 
better performance with p < 0.05 on paired t-test, whereas symbol § refers to a 
statistical significance of p < 0.01 . 



Table III 



Seamentation 
Method 


Lumen 


Plaaue 




Media 


AD (mm) 


HD (mm). 


AD (mm) 


HD (mm) 




AD (mm) 


HD (mm) 


FMM-3D 


0.072 


0.226 


0.061 


0.154 




0.063 


0.164 


PDFs . 


±0.062 


±0.074 


±0.038 


±0.046* 




±0.038 


±0,048§ 


FIVIM-3D 


0.069 


0.197 


0.060 


0.173 




0.063 


0.180 


Gradient 


±0.056 § 


±0.085 § 


±0.044 


±0.050 




+0.044 


±0.052 



[0149] The average and Haussdorf distance were chosen as 

comparison metrics instead of area or perimeter differences because they directly 
depict point-to-point contour variations. As can be seen in Table III. very low 
average and Haussdorf distances values were obtained, for both PDF- and 
gradient-based three-dimensional (3D) fast marching, demonstrating that the 
method is very powerful for simulated B-mode IVUS segmentation. In fact, average 

I 

deviation ranged from 0.060 to 0.072 mm and worst point-to-point distances were 
between 0.154 and 0.226 mm, which Is highly satisfactory. Lower Haussdorf 
distances were obtained on lumen boundary with the gradient method (p < 0.01) 
because the blood and intima Interface generally produces bright echoes for which 
the gradient infbmiatlon is significant. However, on less contrasting boundaries 
such as the intima (plaque) and media interfaces, statistically significant lower 
Haussdorf distances (p < 0.05) were achieved with the PDF-based method. 

[01501 Examples of results obtained with the gray level Rayieigh PDF 

method and with the gray level gradient method for the in-vivo IVUS data are 
displayed In Figures 9a to 10c. The lumen, intima and media detected boundaries 
are presented for a first cross-sectional IVUS image (Figures 9a) and a second 
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different cross-sectional IVUS- image (Figure 10a). 

[0151] In Figures 9b and 10b, the detected boundaries of tiie lumen, 

the intlma and media are estimated with the gray level PDF based fast-marching 
method and, in Figures 9c and. 10c, with the gray level gradient based fast- 
marching method. 

[0152] A qualitative analysis of the gray level PDF and gray level 

gradient fast marching segmentation methods reveals that the detected 
boundaries are very close to all vessel layers. More specifically, Figures 9a to 10c 
show that vessel layer boundaries of in-vivo IVUS images can be Identified even if 
the contrast is very low, as Illustrated at 4 o'clock for the collateral vessel In Figure 
9a. Also, detected boundaries on Figure 10b and 10c demonstrate that non- 
circular lumen may be detected with fast marching methods. 

[01 53] The following Table IV shows the average distance (AD) and the 

Haussdorf distance (HD) between detected boundaries on in-vivo data for the gray 
level PDF and gray level gradient fast marching methods for different manual 
initializations of the interfaces. 



Table IV 



Seamentation 
Method 


Lumen 
AD (mm) HD (mm) 


Plaaue 
AD (mm) HD(mm) 


Mec 
AD (mm) 


Jia 

HD (mm) 


FMM-3D 
PDFs 


0.092 
±0.095 • 


0.270 
±0.147 * 


0.092 
±0.078 


0.256 
±0.102 § 


0.092 
±0.083 


0.256 
±0.113* 


FMM-3D 
Gradient 


0.092 
±0.104 


0.317 
±0.148 


. 0.090 
±0.080 * 


0.287 
±0.092 


0.085 
±0.088 * 


0.302 
±0.107 



[0154] In Table IV, FMM refers to the fast-marching segmentation 

method, symbol * indicates a statistically significant better perfonnance with p < 
0.05 on paired t-test. whereas symbol § means a statistical significance of p < 
0.01. 
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[0155] Results indicate that gray level PDF fast marching has the 

smallest Haussdorf distances (p < 0.01), which remains under 0.270 mm for all 
boundaries' compared to a value of up to 0.317 mm for the gray level gradient 
implementation. PDF fast marching also has relatively small average distances 
between contours, of 0.092 mm and lower, but are significantly higher than Intima 
and media average distances obtained with the gray level gradient method (p < 
0.05). However, the differences between these distances are generally small 
(lower than the pbcel size). Thus, 3D fast marching detected boundaries have small 
variations when initialized differently and the maximum distance to the closest 
point, representing the worst case, generally stays low. This tends to indicate that 
the segmentation perfomiance is good in regions lacking information, for example 
when the boundaries to be detected are covered v^th catheter ring-down artifacts 
of lost behind calcium deposits. 

[0156] Figure 11 shows a 3D reconstruction of the lumen and media 

contours obtained with gray level PDF 3D fast marching segmentation for which a 
double stenosis in that patient is clearly seen. In the figure, the light gray 
corresponds to the inner portion of the media, whereas the darl< gray is the vessel 
lumen. The gray^level gradient fast marching method provided similar qualitative 
results (data not shown). . 

[0157] As mentioned hereinabove, the PDF-based fast-marching 

segmentation method can further exploit RF data in place of B-mode data. On RF 
data, the EM algorithm generally searches for a mixture of Gaussian PDFs 
describing the different layered structures of the vessel wall on IVUS images. 

[0158] Figure 12a shows a simulated 2D RF image taken from the 

whole 3D data set, whereas Figure 12b presents an example of segmentation 
results obtained with the PDF method applied on RF IVUS data. Qualitatively, 
similar perfomiance can be appreciated when com|i)aring those results to Figure 
8c. However, quantitatively, better accuracy is generally expected because of the 
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higher resolution of RF images when compared to B-mode data. 

[0159] In a preliminary version of the IVUS fast marching segmentation 

method disclosed by Boy Cardinal et al,, In "Intravascular ultrasound image 
segmentation: A fast marching method' {Lecture Notes In Computer Sciences. 
Proceedings of MiCCAi 2003: Me6icai image Computing & Computer Assisted 
Intervention, vol. 2879, -2003. pp. 432-439), a 2D version of the fast marching 
method was Implemented. 

[0160] Generally spealdng, a 2D IVUS algorithm uses segmentation 

from previous IVUS images of the catheter pullback to correct initial Interfaces. The 
2D segmentation model disclosed by Roy Cardinal in the above-mentioned study 
was applied to a small IVUS catheter pullback of 200 images. Depending on the 
IVUS application, any dimensions, can be considered for Implementing the fast- 
marching PDF or gradient based method. The present multi-dimensional method is 
general and conceptually considers 1D to ND dimensions, where N Is the order of 
the method. Note that N = 4 considers time varying 3D IVUS data. 

[0161] Since a bigger IVUS B-mode clinical database was available for 

the present study, the 2D version of the fast marching segmentation was applied to 
all available catheter pullbacks. The 2D ImplementaUon of fast marching arrival 
time (from Equation 10) and speed functions (from Equations 11 and 12) is 
generally stralghtfonwartJ. In 2D, 8-connected pixels (26 connected pb(els for the 
hereinabove presented method) were used for averaging neighbors In the speed 
function calculatton. 

[0162] The following Table V shows the average distance (AD) and the 

Haussdorf. distance (HD) between boundaries of the detected vessel layers, from 
different Initializations with 2D fast-marching segmentation. As for the 3D fast 
marching method, the results are for automatically detected PDF- and gradient- 
based algorithms. 
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Table V 



Seamentatlon 
• Method 


Lumen 
AD (mm) HD (mm) 


Plaaue 
AD (mm) HD (mm) 


Media 
AD (mm) HD (mm) 


FMM-2D 
PDFs 


0.093 
±0.096 


0.279 
± 0.149 


0.093 
± 0.078 


0.262 
± 0.102 


0.091 
± 0.083 


0.262 
±0.113 


FMM-2D 
Gradient 


0.096 . 
±0.106 


0.316 
±0.147 


0.095 
± 0.088 


0.299 
±0.100 


0.085 
±0.090 


0.304 
±0.104 



[0163] A two way analysis of variance with pairwlse multiple 

comparisons using the Turkey test was performed on average and Haussdorf 
distances for 2D (Table V) and 3D (Table IV) fast marching. 

[0164] . Statistical tests showed that average distances from the 2D fast 

marching detected boundaries in Table V were not different from the 3D fast 
marching results for all blood vessel Igyers. It should be noted that 2D algorithms 
used segmentation from previous images of the catheter pullback to conrect initial 
contours, which Increased boundary precision. Thus, alternatively, the 
segmentation performance can be Increased by combining this type of conrection 
strategy in the 3D fast marching method, by using a multi-scale segmentation 
approach to Initialize a higher resolution data set with low resolution segmentation 
results of the same catheter pullback. As for the 3D fast marching method, the- 
gray level PDF fast marching in 2D . had lower Haussdorf distances tiian the gray 
level gradient method (p < 0.05). Since good average distance perfomiance was 
achieved with the gray level gradient method in both 2D and 3D fast marching, this 
information can be added with advantage to the PDF speed function of Equation 
11. 

[0165] A second non-restrictive illustrative embodiment of the method 

and device according to the present invention will now be described. For tiie sake 
of brevity, only the differences between the method and device according to the 
first non-restrictive illustrative embodiment and the second non-restrictive 
illustrative embodiment are described hereinbelow. 
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[0166] In the secxjnd non-restrictive illustrative embodiment, the fast 

marching niethod has been modified to enhance the efficiency In computing time. 
In the fast marching method In accordance with the first Illustrative embodiment, 
several Interfaces simultaneously propagate across the IVUS Image data. In the 
propagation process, the' propagating interfaces and their neighboring areas are 
explored in a significantly detailed manner. Generally, in the segmentation method 
as described In the first illustrative embodiment, all pixels are analyzed and the 
propagation process takes into account all preceding interface neighbors through 
the arrival time map construction. Computation time Is thus increased as the initial 
interface propagates in a larger Initial segmentation area. 

10167] In ttie first Illustrative embodiment, ttie position of the initial 

interfaces is calculated from a shrunk version of manually or automatic traced 
initialization contours taken along a longitudinal plane. Examples of tiie position of 
initial interfaces calculated from initialization contours are shown in Rgures 13a 
and 13b. ' • 

[0168] The black region represents the unexplored propagating area 

130, the gray pixels on each side of the propagating area 130 con-espond to the 
propagating interfaces 132,134 and the anx)ws 136a,1 36b. 136c and 136d 
represent the propagation direction of ttie propagating interfaces 132,134. The 
dashed line 138 represents the desired boundary to reach and the solid line is an 
example of initial interface 140 from which tiie initial propagating interfaces 132, 
134 were calculated. In Figure 13a, tiie interfaces will not detect the boundary 138 
because the propagation area 130 was not set wide enough to completely include 
tiie boundary 138, as in Figure 13b. However, propagation in Figure 13a will be 
completed faster ttian tiie case of Rgure 1 3b. 



[0169] • Therefore, to decrease computational load, shrinking can be 
diminished to create a smaller propagating area 130 (Figure 13a). However, 
because tiie fast marching method propagates an interface 132,134 under a 
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unidirectional speed function (see anx}ws 136a,136b or 136c,136d) the boundary 
138 to be detected must be iocated inside the propagating area 130 that wiil be 
explored during the propagation. 

[0170] A compromise between the dimension of the propagating area 

130 and the computation time is sought. With loiown 2D fast marching 
segmentation method, this problem is generally solved by using segmentation 
results from previous 2D images of the catheter pullbacic to con-ect initial 
Interfaces. The initial Interfaces 140 are then more precise and the propagating 
area 130 can be set with smaller dimensions. 

[0171] in the 3D fast marching segmentation method of the first non- 

restrictive illustrative embodimient, a cbn-ection similar to this 2D correction 
principle can be made through a multi-resolution or multi-scale representation and 
segmentation of the IVUS data. An example of such multi-resolution images in 
IVUS data is shown in Figures 14a-14d, where Mower resolution images are 
obtained by undersampiing flie original IVUS Image by 2' (Operation 191 of Rgure 
20) where /= 3, 2, 1, 0 is tiie resolution, levd con-esponding to Fig. 14a, 14b, 14c 
and 14d, respectively. High scale stiiictures, such as for example the lumen, are 
generally emphasized on lower resolution images. 

[0172] The segmentation results" of a lower resolution representation of 

ttie IVUS data are mapped into the next level of resolution (Operation 192 of 
Figure 20). Boul<em)ul et al., In "Segmentation of ulti^sound Images - 
multiresolutlon 2D and 3D algorithm based on global and local stati'stics" (Pattem 
Recognition Letters, 24:779-790, 2003) and Mignotte et al„ In "A multiscale 
optimization approach for the dynamic contour-bksed boundary detection issue" 
(Computerized Medical Imaging and Graphics, 25(3):265-275, 2001) propose 
related concepts. 



[0173] 



These segmentation-mapped results are used to initialize the 
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interface propagation at lliis higher resolution level (Operation 193 of Figure 20). 
At a low-resolution level, a fast coarse exploration of a wide propagating area is 
performed to bring the propagating interfaces 132,134 closer to the desired 
boundaries 138 (Operation 201 of Figure 21). The propagation area 130 can then 
be reduced at each higher resolution level since the Interfaces 132,134 are 
iteratlvely corrected (Operation 202 of Figure 21). A larger propagation area 130 
can thus be explored in less calculation over all possible resolutions. 

[0174] At a resojutipn level /, a pixel represents a 2' x 2' block of 

pixiels from the driglnal resolution image. In Figures 14a to 14d, the undersampling 
that was used resulted In loss of InformaBon at low resolution levels. To overcome 
this loss of infomiatlon, a multlscale PDF-based velocity function (similar to 
Equation 11) Is developed and used (Operation 203 of Figure 21), where P(l8) is 
replaced with the lil<elihood of the 2' x 2' block of pixels coresponding to Is when 
the propagation Is done at level /. as given by equation 13 : 



[01751 P(7,)=nP(/„) (13) 



[0176] where bi is the block of 2' x 2 ' pixels and P( ; Is the occurring 

probability of the gray level value of pbcel s/ In the zero resolution Image / . 

[0177] The multiresolution and multiscale fast marching segmentation 

methods of Figures 20 and 21 generally allows to iteratlvely Improve the accuracy 
of the detected boundaries without increasing the computation time. 

[0178] A third non-restrictive illustrative emjjodlment of the method and 

device according to the present invention will now be described. For the sake of 
brevity, only the differences between the third; non-restrictive Illustrative 
embodiment and the first non-restrictive illustrative erjibodiment will be described 
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hereinbelow. 

[0179] In this" third non-restrictive IllustratVe embodiment, the fast 

marching method has been modified to automatically find the initial interlaces for 
the layers of the vessel wall (lumen, inside and outside contours of the media) by 
using likelihood maps of. each components of the vessel wall (lumen, intima and 
plaque, media, and sun^ounding tissues), which are calculated according to the 
estimated PDF mixture parameters. This approach can be. seen as an altematlve 
to the manual initializations of the vessel interfaces described hereinabove. 

[01801 "ThQ initialization procedure generally finds a rough estimate of 

the true boundaries of the layers that will be further refined into accurate wall 
contours with the multiresoiuBon or multiscale gray levej PDF fast marching 
method. . . ' 

[0181] The initialization procedure generally starts on a small subset of 

contiguous 2D IVUS frames from the whole catheter pullback (Operation 211 of 
Figure 22). For instance,- the subset may contain 2D INAJS frames to get as 
much infonnaaon as possible {Nm = ten (10) for example), while keeping enough 
con-elation between the 2D IVUS flames. For better ri9sults, the selected 2D IVUS 
frames are generally of good quality, with no. calcification shadows and with a 
generally homogeneous plaque, in order to maximize the available Infonmation for 
detennining the initial interface. 

[0182] To find these good quality 2D IVUS frames, a degree of.fitting is 

first calculated between each of the Individual frame histogram and the pullback 
PDF mixture (see Figure 3). The degree of fitting can be measured by the 
Kolmogorov-Smimov test that calculates the distance between the PDF mixture 
and the normalized histogram of an IVUS 2D frame. This test is used to determine 
if an empirical distribution (IVUS histogram) differs from a theoretical distribution 
(mixture). The Kolmogorov-Smimov distance, between the subset init of Nm 
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images and the global IVUS data PDF P. that should be minimized is : 



[0183] i:=ZZl^-Wl 04) 



[0184] where U is the gray level value of pixel s In Image /; kjCQ Is the 

number of pixels having the gray level value U in they* image of the subset; and N 
is the number of pixels In image /. 

[0185] The contiguous frames having the smallest Kolmogorov- 

Smlmov distances are generally chosen. Since the mixture, is calculated over the 
whole catheter, pullback, it represents the average lumen and blood vessel wall. 
This test thus generally selects the frames that are similar to the average catheter 
pullback, and these frames are used to start the calculation of the Initial interface. 

[0186] Calculation of the Initial Interface is Initiated with the search of 

an inner generally elllptlcaj ring shaped region corresponding to the media 
structure 150 (Operation 212 of Figure 22). 9S shown In Figure 15. This media 
region 150 is distributed according to the hereinabove disclosed media PDF and 
. enclosed between a surrounding tissue region 152 and a plaque region 154. The 
media region 150. the surrounding tissue region 152 and the plaque region 154 
are generally of fixed size to simplify the InltiaRzatlon procedure and because only 
a rough estimate of the media is generally necessary. The surrounding tissue 
region 152 and the plaque region 154 do not have to represent the whole plaque 
and surrounding tissues since they are defined to provide additional Infonnatloh 
about the media region 150. 

[0187] The Initialization procedurei generally beglns.wlth the search of 

the media region 150 of the blood vessel because .it Is believed that the elliptical 
constraints. are more easily assumed for tWs layer. It was indeed reported by 
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Hagenaars et al., in "Gamma radiation induces positive vascular remodeling after 
balloon angioplasty: A prospective, randomized intravascular ultrasound scan 
study" (Journal of Vascular Surgery, 36(2):31 8-324, 2002) In which 15 patients out 
of 16 had a dissection after angioplasty of the femoropopliteal artery making the 
lumen in^egulariy shaped. 

[0188] Since IVUS data acquisition is often conducted in 

atherosclerosis treatment trial in which patients undergo angioplasty, the 
in-eguiarities of the layers should be taken into account. Also, the search of the 
initial lumen interface is reduced to the inside region of the media (Operation 213 
of Figure 22), which generally prevent the propagating interface from leaking into 
collateral branches when they are present. Moreover, the elliptical shape of the 
initial media region 150 generally produces a closed initial interface, even if the 
media Is hidden behind a calcification shadow. 

[0189] In order to find the media region 150 in the subset of initial IVUS 

frames, an energy function must be associated with the template 158 of Figure 15. 
In the defomiable model framework disclosed by Jain et al., in "Defomnable 
template models: A review" (Signal Processing, 71(2):109-129, 1998) and by 
Zhong et al., in "Object tracking using deformable templates" (In Sixth Intemational 
Conference on Computer Vision; pages 410-445, 1998), the energy function is 
generally defined. to be minimal when the template fits the searched region. 
Defonnations are applied to the template 158 to achieve a minimum energy 
function on the image. 

[0190] The energy fundlon that should ^)e mlnimlzal to find the media 

region 1 50 is given by the following Equation 1 5: 



[0191] 



1 Nttiit 
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[0192] where P(ll |r/) is the occurring probability of the gray Jevel value 

of pixel s in the J ^ IVUS image of the IniHal' subset according to the PDF of region 
n; R^lr^yr^yr^] are the plaque region (r^) 154, media region (r^) 150 and tissue 

region (r^) 152 of the template 158; Ninu is the size of the initial frame subset. 

10193] Generally only linear transformations such as for example 

translations, stretchings^ and rotations are applied to the template 158 because 
only a generally rough estimate of the media region 150 is needed. The template 
158 fitting Is perfomied at a reduced resolution level / = 1, as described in the 
secopd non-restrictive illustrative embodiment, In order to minimize the 
computation time, while keeping a large enough media to wori^ with. Different 
known minimization algorithm can be used for minimizing the defomiatlon model. 

[0194] The lumen region (not shown in Rgure 15) Is then searched 

from the defined initial media boundary. A region without geometrical shape 
constraints is computed: the lumen likelihood map is calculated and a lumen 
region is grown or propagated using this map (Operation 214- of Figure 22). The 
lumen region generally starts at the IVUS catheter that is inside the lumen and 
generally located in the vicinity of the center of all puilback frames. 

[0195] The lumen region grows by adding the pixels that are most likely 

to be inside the liimen according to the occurring probability, for example if the log- 
likelihood -fllogP(// 1 r/) of pixel s in the Nmn image subset according to the PDF 
>> 

of the lumen region n is low enough. The. region is generally fort)idden to grow 
beyond the boundary of the media region 1 50. 

[0196] The media and lumen regions are then adjusted or fitted to the 

Nfnit contiguous 2D IVUS frames that were used Jn the Initial subset (Operation 215 
of Figure 22). Linear transfomiatlons may be applied to these media and lumen 
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regions to maximize their likelihood to each of the contiguous 2D IVUS frames. 
This step is usually performed as for the herelnat)ove described media defomnable 
template, but ttie Initial intlma and plaque region 154 is . bounded by the lumen 
region and the media region 150, and the initial surrounding tissues region 152 is 
bounded on only one side by the media region 150. 

[0197] This procedure is generally repeated for the next subset of 

contiguous 2D IVUS frames in the catheter pullback. However, the process for 
each contiguous 2D IVUS frame generally starts with the results of the previous 
defined media template 156. The growth of the lumen region generally starts from 
a shrunk version of the previous average lumen region. The whole IVUS image 
pullback is therefore initialized in that manner. 

[01981 Alternatively, the segmentation fast marching method as 

described in ttie first illustrative embodiment may use a combination of the gray 
level gradient Information and the gray level PDF Information in the calculation of 
the initial Interfaces if using only the gray level PDF information turns out to be 
Insufficient to generate an automatic Initialization as described In tiie third 
illustrative embodiment. The gray level gradient Infomiaflon could also be 
integrated to the Interface velodty function of Equation 1 1 . 

[Oi^S] In the case of very low-quality images or high ultrasound 

attenuation limiting penetration within flie vascular wall, ttie proposed initial 
boundary calculation procedure might fail to find some Initial contours or region 
boundaries. For these particular cases, minimal user interaction might be required 
to con-ect some regions of tiie Interfaces. If necessary, this interaction may further 
be included In flie segmentation process to minimize tiie occurrences of having to 
re-segment accurately found boundaries. 

[0200] In the case where a single boundary is available, such as for 

example the boundary between the lumen and the intima, the elliptical template 
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may be modified In such a way as to provide a two (2) region template generally 
corresponding to the regions between the single boundary. The energy function of 
equation 15 thus becomes generally restrained to two (2) regions and the 
remaining, initialization procedure remains generally similar to the case of the 
multiple boundary initialization. 

/ 

[0201] A fourth, non-restrictive illustrative embodiment of the method 

and device according to the present invention, will now be described. For the sake 
of brevity, only the differences between this fourth illustrative embodiment and the 
first illustrative embodiments will be described. 



[0202] In this fourth non-restrictive illustrative embodiment, the fast 

marching method has been modified to replace the EM local algorithm presented 
in the first illustrative embodiment. The ElVI algorithm Is a local algorithm in which 
the- neighbors infomiatlon is missing. This infomiation is generally required, such 
as for example, in the case of heterogeneous plaque where the PDFs are 
generally more difficult to estimate. In addition, convergence is generally very slow 
with the EiVI algorithm such that It can take a significant number of iterations in 
order to be able to estimate the mixture parameter © of Equation (2). 

[0203] The automatic initial contour detection procedure presented in 

the previously presented thlitl illustrative is based on the IVUS PDF infomnation. In 
cases where the EM .algorith.m cannot be used, the iterative conditional estimation 
(ICE) algorithm that was previously proposed for the mbrture parameter estimation 
of incomplete data by Pleczynski In "Champs de mari«ov caches et estimation 
condiaonnelle iterative" (Traltement du Signal. 11 (2):141 -153. 1994) generally 
represents a more robust algorithm, which generally converges faster than the EM 
algorithm. • 



[0204] In the PDF mixture estimatfon presented in the first illustrative 

embodiment, the random variables (X, Y) are refen-ed to as the complete data 
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where Y is the gray level taking values In [1, 256] (observed data), and X Is the- 
tissue label taking values [1, M] for a PDF mixture modeling M different tissues 
(hidden information). For the set of pixels S, the realization y-{y^X^ of Y Is the 

INAJS B-mode or RF Image and x=^{x^\^ are the unknown pixel labels. The EM 

algorithm is considered local because the labels are considered independent. In 

the ICE algorithm. X is suppoised l\/1aricovlan I.e. Px (x) Is defined with respect to 
the following neighborhood energy function: 



[0205] 



P(x)=exp 



(16) 



[0206] where ^ is an energy function and the summation Is for all pairs 

of pixel neighbors (sj). 

[0207] The first operation of the ICE algorithm is to simulate, such as 

for example with the Gibbs sampler, n realizations (ji:i,...,;c„) of X according to the 

posterior distribution iA'|y,0(^>',0')» with ©' the initial or the previous iteration 

estimate of the PDF mixture parameter 0 (Operation 221 of Figure 23). The 
posterior distribution is computed using Bayes.rule from the known Px(x) and 
Py|x;6>fy|x> ®'); and the neighbortiood Vs of pixel s : 



[0208] 



(17) 



[0209] 



With these simulations of the hidden data, n sets 
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{ixi,y),...Xx„,yj) of complete data are available. 

[0210] The next operation (Operation 222 of Figure 23) is to calculate 

the new value of ®' =-\^{xt,y)+...+®{x„,y)], where 0 is a parameter estimator 
n 

of the complete data (maximum likelihood for example). Operations 221 and 222 of 
Figure 23 are generally repeated until convergence of the mixture parameter 
estimate is achieved (Operation 223 of Rgure 23): 

[0211] For the Rayleigh mixture, it is assumed that each layer structure 

of the B-mode images is a generally unifomn scattering tissue with a significantly 
large number of diffusers because the Rayleigh PDFs model the gray level 
distribution of the ultrasound signal under that condition, as disclosed by VVagner 
et al., in "Statistics of speckle in ultrasound B-scans" (IEEE Transactions on.Sonics 
and Ultrasonics, 30(3): 156-1 63. 1983). The same reasoning applies to Gaussian 
PDFs describing RF IVUS Images. 

[0212] In the case of highly heterogeneous plaque layer of a diseased 

patient, the Rayleigh or Gaussian PDF might not be sufficient to model the pixel 
gray level distribution. Distributions other than Rayleigh or Gaussian have been 
Investigated. In modeling of the ultrasound B-mode envelopes or RF signals, 
respectively: RIcian .distribution as disclosed by Wear et al., In "Statistical 
properties of estimates of signal-to-nolse ratio and number of scatterers per 
resolution cell" (Journal of the Acoustical Society of America. 102(1):635641. 
1997), K distribution as disclosed by Dutt et al.. in "Statistics of the log-compressed 
echo envelope" (Journal of the Acoustical Society of America, 99(6):381 7-3825; 
1996) and NakagamI distribution as disclosed by Shankar In "A general statistical 
model for ultrasonic backscattering from tissues" (IEEE Transactions on 
Ultrasonics. Fem)electrics, and Frequency Control, 47(3):727- 736. 2000). 



[0213] 



The ICE algorithm generally has no limitation for the type of 
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statistical distribution to be modeled, as long as a parameter estimate of the 
complete data can be calculated. Moreover, if a model of mixed distribution types 
is necessary, the generalized ICE algorithm (GICE) as disclosed by-Delingnon et 
al., In "Estimation of generalized mixtures and Its application In image 
segmentation" (IEEE Transactions on Image Processing, 6(10):1364-1375, 1997) 
can be used. i3ICE generaliy.provides parameter estimates for mixtures composed 
of a various number and type of statistical distributions. 

[0214] A fifth non-restrictive illustrative embodiment of the method and 

device according to the present invention will now be described. For the sake of 
brevity, only the differences between the fifth and first illustrative embodiments will 
be described hereinbelow. 

[0215] In this fifth non-restrictfve illustrative embodiment, the 

segmentation fast marching method allows to treat and analyze, In addition to the 
volumic Infomiatlon obtained from the boundary layer detections of the blood 
vessels, dynamic data retrieved from IVUS pullbacks defining a fourth dimension. 
The dynamic data generally relates to the cyclic pulsation occurring in the blood 
vessels. 

[0216] The cyclic variations of the vessel dimensions combined to' 

cardiac motion (for coronary IVUS) was described, in the literature, as the 
sawtooth artifact which is generally visible on longitudinal view of the IVUS volume 
and generally caused by the blood vessel pulsations. As shown by the arrow in 
Rgure 16, the sawtooth artifact Is also present In IVUS data of fenroral arteries 
even without cardiac motion. In common femoral arteries for Instance, diameter 
measurements generally , vary from 6.8 inm in diastole to 7.2 mm in systole for 
patients with lower limb peripheral vascular disease, as disclosed by Tal et al. In 
"In vivo femoropopliteal arterial wall compliance In subjects with and without lower 
limb vascular disease." (Journal of Vascular Surgery, 30(5):936-945, 1999). 
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p)217] Electrocardiogram-gating (ECG-gating) acquisition was 

proposed by von BirQeien et al., in "ECG-Qated three-dimensional intravascular 
ultrasound, feasibility and reproducibility of the automated analysis of coronary 
lumen and atherosclerotic plaque dimensions in humans" (CircuIatlon,-96(9):2944- 
2952, 1997) to remove this artifact. This is generally accomplished by acquiring 2D 
IVUS frames at a precise moment of the cardiac cycle, commonly at the end of 
diastole, which generally gives more accurate and reproducible volumic 
measurements, as disclosed by von Birgelen et al. and by Bruining et al., in "ECG- 
gated versus nongated three-dimensional intracoronary ultrasound analysis:, 
implications for volumetric measurements" (Catheterization, and Cardiovascular 
Diagnosis. 43:254-260. 1998). 

[0218] Because ECG-gating hardware Is generally not widespread, 

retrospective gating was proposed to remove the cydjc changes on non-gated 
IVUS puliback. Change tracking In semi-automatically detected iumen contour was 
first proposed by Nadkamni et a!., in "Image-based retrospective cardiac gating for 
three-dimensional intravascular ultrasound Imaging" (SPIE Proceedings: Medical 
Imaging: Ultrasonic Imaging and Signal Processing, volume 4687, pages 276- 
284, 2002). 

[0219] Another method searched for cyclic variations in contour 

features calculated In a pre-processing step as disclosed by de Winter et al., in 
"Retrospective Image-based gating of Intracoronary ultrasound images for 
Improved quantitative analysis: The intelligate method" (Catheterization and 
Cardiovascular Diagnosis. 61:84-94, 2004). The most recent retrospective gating 
proposed method Is based on variations of the Images mean gray level values by 
ZhU et .al., in "Retrieval of cardiac phase from IVUS sequences" (SPIE 
Proceedings: Medical Imaging: Ultrasonic Imaging and Signal Processing, volume 
5035. pages '135^146. 2003) which states that the bigger systolic lumen, that is 
hypbecholc, generally decreases the mean gray level value of the image. 
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[0220] Some measurements can be made from the cyclic vessel 

variations. It was demonstrated by Shaw et al., In "Detemnlnants of coronary artery 
compliance in subjects with and without angiographic coronary artery disease" 
(Joumal of the American College of Cardiology, 39(1 0):1 637-1 643, 2002.) that 
plaque compression is related to the vessel cross-sectional compliance. Also the 
lumen cross-sectional area (CSA) difference between systolic and diastolic 
measurements was significantly greater in yellow plaque which generally consists 
of thin, fibrous cap with Jipid-rich core and inadequate collagen content), than white 
plaque which consists of thick fibrous cap or completely fibrous, as disclosed by 
Takano et aL, in ''Mechanical and structural characteristics of vulnerable plaques: 
Analysis by coronary angioscopy and intravascular ultrasound" (Journal of the 
American College of Cardiology, 38(1 ):99-104, 2001 ). * 

[0221] Thus, the cyclic pulsation contains information about volumic 

changes of the blood vessel wall that is generally lost when the acquisition is ECG- 
gated. The vessel pulsation infomriation from non-gated acquisition may be kept 
and used to reconstruct the Vessel wall in 3D, at different moments of the cardiac 
cycle. With this fourth-dimensional reconstruction of the vessel wall, volumic 
accuracy and reproducibility can be achieved for measure;nents made on 3D 
image sets at specific moments of the cardiac cyde. 

[0222] To perform 4D reconstruction of the blood vessel wall, detected 

boundaries from each 2D IVUS frames first have to be classified in different wall 
pulsation phases. This step may be achieved by searching periodic components In 
measurements calculated from the detected boundaries. 

[0223] Figure 17 shows the lumen area calculated from the 

segmentation results obtained from the semi-automatic l='DF-based fast marching 
method as described hereinabove. A cyclic variation is visible and could be 
assessed on these retrospective measurennents. However, it Is also possible to 
include the wall pulsation assessment in the boundary detection process, such that 
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this information is used to heip the segmentation of the IVUS pullback. 

[0224] The wall pulsation assessment may be initiated during the initial 

contour calculation procedure (Operation 231 of Figure 24), and refined when the 
boundary detection Is finished. The wall pulsation Is then generally divided In a 
discrete number of phases (Operation 232 of Figure 24); such as for example 
^stole, beginning of diastole and end of diastole, and a label Is assigned to each 
wall pulsation phase (Operation 233 of Figure 24). The pullback 2D IVUS frames 
are tiien classified and assigned to Oie corresponding pulsation phase label 
(Operation 234 of Figure 24). 

[0225] As shown In Figure 17, Oie lumen area contains the cydlc 

pulsation Information. Further, the area variation between adjacent frames is 
generally used to define the classification. Since tiie pulsation Is usually periodic, 
the cyclic pulsation infonnation is used to defomi the initial subset template 158 
according to variations of the expected pulsation in tiie Initial contour calculation of 
individual IVUS images (Operation 215 of Rgure 22). The assignment of the 
pulsation phase labels to following frames can also take advantage of the periodic 
infonnation. 

[0226] At the end of the iniflaiization process, each 2D IVUS frame is 

Identified with a wall pulsation phase label. However, tiiese labels may change 
because, at the end of the segmentation process, more accurate lumen areas are 
calculated. The Initial labels iare therefore adjustable according to their initial value, 
to "the variations in area difference measurements and to Hie expected value 
according to the periodic variation (Operation 235 of Figure 24). 

[0227] With tills pulsation assessment, tiie 4D data set are divided In 

3D data sets composed of all IVUS 2D frames associated to a spedfip cardiac 
phase label and corresponding to ttie different phases of tiie blood .vessel's 
pulsation (Operation 236 of Rgure 24). Volumic (boundary) measurements can 
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then be made on tfiese separate 3D data sets (Operation 237 of Rgure 24). 

[0228] Although the present Invention has been described hereinabove 

by way of non-restrictive illustrative embodiments thereof, It can be modified at will, 
within the scope of the appended claims without departing from the spirit and 
nature of the subject invention. 
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